/*     CalculiX - A 3-dimensional finite element program                 */
/*              Copyright (C) 1998-2021 Guido Dhondt                          */

/*     This program is free software; you can redistribute it and/or     */
/*     modify it under the terms of the GNU General Public License as    */
/*     published by the Free Software Foundation(version 2);    */
/*                                                                       */

/*     This program is distributed in the hope that it will be useful,   */
/*     but WITHOUT ANY WARRANTY; without even the implied warranty of    */ 
/*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the      */
/*     GNU General Public License for more details.                      */

/*     You should have received a copy of the GNU General Public License */
/*     along with this program; if not, write to the Free Software       */
/*     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.         */

#ifdef __WIN32
_set_output_format(_TWO_DIGIT_EXPONENT);
#endif

#ifdef CALCULIX_MPI
#include <spoolesMPI.h>
#endif

#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include "CalculiX.h"

#ifdef CALCULIX_MPI
ITG myid = 0,nproc = 0;
#endif

struct timespec totalCalculixTimeStart,totalCalculixTimeEnd; 

void CalculiXstep(int argc,char argv[][133],ITG **nelemloadp,double **xloadp,
		  ITG *nload,char **sideloadp,double *timepar,ITG *ne,
                  ITG **ipkonp,ITG **konp,char **lakonp,ITG *nk,double **cop,
                  double **voldp,double **veoldp,double **accoldp,
                  ITG *nboun,ITG **ikbounp,ITG **ilbounp,double **xbounp,
                  ITG *nmethod,char *externalfaces,double *delexternalwork,
                  ITG *inputsteps,ITG *iperturb,ITG *irstrt,char **filabp,
                  ITG *nlabel)
{

  /* in FORTRAN convention:

     nelemload(1,i) (ITG): element for loading i
     nelemload(2,i) (ITG): node for the environmental temperature 
                           (heat transfer analysis)

     xload(1,i) (double): magnitude of load at end of a step or, for heat 
                          transfer analyses, the convection or the radiation 
                          coefficient
     xload(2,i) (double): the environment temperature (only for heat transfer 
                          analyses)
     
     nload (ITG): number of distributed loads

     sideload(i): load label; indicates element side to which load is applied

     timepar(1) (double): time increment
     timepar(2) (double): step length
     timepar(3) (double): minimum time increment allowed
     timepar(4) (double): maximum time increment allowed
     timepar(5) (double): time increment for CFD-calculations

     ne (ITG): highest element number
     ipkon(i) (ITG): pointer for element i into field kon
     kon(ipkon(i)+1,....ipkon(i)+nope) (ITG): topology (node numbers) for 
                                              element i; "nope" depends on the 
                                              element type (stored in lakon)
     lakon(i) (character*8): label for element i

     nk (ITG): highest node number
     co(1..3,i) (double): coordinates of node i
 
     vold(0,i) (double): temperature in node i
     vold(1,i) (double): displacement in global x in node i
     vold(2,i) (double): displacement in global y in node i
     vold(3,i) (double): displacement in global z in node i

     veold(0..3,i) (double): first time derivative of field vold
     accold(0..3,i) (double): second time derivative of field vold

     nboun (ITG): number of single point constraints (SPC)

     ikboun(i) (ITG): sorted degrees of freedom corresponding to the SPC's. 
                      The degree of freedom corresponding to global direction 
                      k in node j is idof=8*(j-1)+k. Applicable directions 
                      are 0 (temperature), 1 (global x-dir), 2 (global y-dir) 
                      and 3 (global z-dir)

     ilboun(i) (ITG): SPC number corresponding to the sorted degrees of freedom
                      in ikboun

     xboun(i) (double): value of the SPC

     nmethod (ITG): nmethod=1: static analysis   
     nmethod=2: frequency analysis  
     nmethod=3: buckling analysis 
     nmethod=4: (linear or nonlinear) dynamic analysis 
     nmethod=5: steady state dynamics analysis 
     nmethod=6: Coriolis frequency calculation 
     nmethod=7: flutter frequency calculation 
     nmethod=8:  magnetostatics 
     nmethod=9:  magnetodynamics 
     nmethod=10: electromagnetic eigenvalue problems 
     nmethod=11: superelement creation or Green function calculation 
     nmethod=12: sensitivity analysis  
                   
     externalfaces: (81 characters) surface name of the faces in contact with 
                    the fluid; for this surface the work will be calculated at 
                    the end of the present routine before returning the control
                    to the calling program. The surface must be defined in the
                    CalculiX input deck as a facial surface, i.e. consisting of
                    faces (not nodes). In the calling program of CalculiXstep 
                    there must be a blank behind this name.

     delexternalwork (double): total external work on the faces contained in 
                               the surface defined by externalfaces performed 
                               during this call

     inputsteps (ITG): number of steps to be read from the input deck. This 
                       requires (inputsteps) calls to CalculiXstep.c. Call 
                       (inputsteps+1) is the first call for which the input 
                       deck is not further read, i.e. in this call the 
                       information from step (inputsteps) is taken, possibly
                       modified by the information transferred by the user into
                       CalculiXstep.c through the argument list.

     iperturb(1) (ITG): perturbation parameter: specifies which nonlinear 
                        procedure is to be selected. Should not be changed by 
                        the calling program
     iperturb(2) (ITG): if 0: geometrically linear calculation (linearized 
                              strain)
                        if 1: geometrically nonlinear calculation (quadratic 
                              terms in the strain are taken into account)
                              can be changed by the calling program

     irstrt(1) (ITG): if   0: no restart (neither READ nor WRITE)
                      if > 0: every irstrt steps the results are stored in file 
                           jobname.rout
                   Remark: reading a restart file is obtained by inserting 
                           *restart,read as first executable statement in a 
                           CalculiX input deck
                 
     irstrt(2) (ITG): if 0: no OVERLAY while writing a restart file
                      if 1: OVERLAY while writing a restart file
                
     filab(87*nlabel): character field telling which fields have to be stored 
                       for which sets (cf. CalculiX manual for more details). 
                       Setting this field to blank suppresses any
                       frd-output

     nlabel: total number of labels.


  */ 


  // start change DLR

  /* the contents of all variables has to be kept between two calls of
     CalculiXstep, therefore they are declared as static */
  
  static FILE *f1;
    
  static char *set=NULL,*matname=NULL,*orname=NULL,*amname=NULL,
    *labmpc=NULL,*prlab=NULL,*prset=NULL,
    jobnamec[792]="",jobnamef[132]="",output[5],*typeboun=NULL,
    *inpc=NULL,*tieset=NULL,*cbody=NULL,fneig[132],*sideloadtemp=NULL,
    kind1[2],kind2[2],*heading=NULL,*objectset=NULL;
  
  static ITG *nodeboun=NULL,*ndirboun=NULL,*ipompc=NULL,
    *nodempc=NULL,*nodeforc=NULL,*ndirforc=NULL,
    im,*inodesd=NULL,nload1,*idefforc=NULL,
    *nactdof=NULL,*icol=NULL,*ics=NULL,itempuser[3],
    *jq=NULL,*mast1=NULL,*irow=NULL,*rig=NULL,*idefbody=NULL,
    *ikmpc=NULL,*ilmpc=NULL,
    *nreorder=NULL,*ipointer=NULL,*idefload=NULL,
    *istartset=NULL,*iendset=NULL,*ialset=NULL,*ielmat=NULL,
    *ielorien=NULL,*nrhcon=NULL,*nodebounold=NULL,*ndirbounold=NULL,
    *nelcon=NULL,*nalcon=NULL,*iamforc=NULL,*iamload=NULL,
    *iamt1=NULL,*namta=NULL,*iamboun=NULL,
    *nplicon=NULL,*nplkcon=NULL,*inotr=NULL,*iponor=NULL,*knor=NULL,
    *ikforc=NULL,*ilforc=NULL,*iponoel=NULL,*inoel=NULL,*nshcon=NULL,
    *ncocon=NULL,*ibody=NULL,*ielprop=NULL,*islavsurf=NULL,
    *ipoinpc=NULL,mt,nxstate,nload0,iload,*iuel=NULL,*ne2boun=NULL,
    *irandomtype=NULL,irobustdesign[3],*iparentel=NULL,ifreebody,
    *ipobody=NULL,inewton=0,*iprfn=NULL,*konrfn=NULL;
     
  static ITG nmpc,nforc,nprint,nset,nalset,nentries=18,
    neq[3],i,mpcfree,mei[4],j,nzl,nam,nbounold,
    nforcold,nloadold,nbody,nbody_,nbodyold,network,nheading_,
    k,nzs[3],nmpc_,nload_,nforc_,istep,istat,nboun_,nintpoint,
    nmat,ntmat_,norien,ithermal[2]={0,0},nmpcold,
    iprestr,kode,isolver,nslavs,nkon_,ne0,nkon0,mortar,
    jout[2],nkon,idrct,jmax[2],iexpl,nevtot,ifacecount,
    iplas,npmat_,mi[3],ntrans,mpcend,namtot_,iumat,iheading,
    icascade,maxlenmpc,mpcinfo[4],ne1d,ne2d,infree[4],
    callfrommain,nflow,jin=0,nener,jrstrt,nenerold,
    nline,*ipoinp=NULL,*inp=NULL,ntie,ntie_,mcs,nprop_,
    nprop,itpamp,iviewfile,nkold,nevdamp_,npt_,cyclicsymmetry,
    nmethodl,iaxial,inext,icontact,nobject,nobject_,iit,
    nzsprevstep[3],memmpcref_,mpcfreeref,maxlenmpcref,*nodempcref=NULL,
    *ikmpcref=NULL,isens,namtot,nstam,ndamp,nef,inp_size,
    *ipoinp_sav=NULL,*inp_sav=NULL,irefineloop=0,icoordinate=0,
    *nodedesi=NULL,ndesi=0,nobjectstart=0;

  static ITG *meminset=NULL,*rmeminset=NULL;

  static ITG nzs_,nk_,ne_,nset_=0,nalset_,nmat_,norien_,nam_,
    ntrans_,ncs_,nstate_,ncmat_,memmpc_,nprint_,nuel_=0;
    
  static double *coefmpc=NULL,*xforc=NULL,*clearini=NULL,
    *xbounold=NULL,*xforcold=NULL,*randomval=NULL,
    *sti=NULL,*xloadold=NULL,*xnor=NULL,
    *reorder=NULL,*dcs=NULL,*thickn=NULL,*thicke=NULL,*offset=NULL,
    *elcon=NULL,*rhcon=NULL,*alcon=NULL,*alzero=NULL,*t0=NULL,*t1=NULL,
    *prestr=NULL,*orab=NULL,*amta=NULL,
    *t1old=NULL,*eme=NULL,*plicon=NULL,*pslavsurf=NULL,*plkcon=NULL,
    *xstate=NULL,*trab=NULL,*ener=NULL,*shcon=NULL,*cocon=NULL,
    *cs=NULL,*tietol=NULL,*fmpc=NULL,*prop=NULL,*t0g=NULL,*t1g=NULL,
    *xbody=NULL,*xbodyold=NULL,*coefmpcref=NULL,*dacon=NULL,*vel=NULL,
    *velo=NULL,*veloo=NULL,energy[5],*ratiorfn=NULL,*dgdxglob=NULL,
    *g0=NULL,*xdesi=NULL;
    
  static double ctrl[57];

  static double fei[3],*xmodal=NULL,alpha[2],ttime,qaold[2],physcon[14];

  static double totalCalculixTime;

  // end change DLR

  
#ifdef CALCULIX_MPI
  MPI_Init(&argc, &argv) ;
  MPI_Comm_rank(MPI_COMM_WORLD, &myid) ;
  MPI_Comm_size(MPI_COMM_WORLD, &nproc) ;
#endif


  // start change DLR

  char *sideload=NULL,*lakon=NULL,*filab=NULL;
  ITG *nelemload=NULL,*ikboun=NULL,*ilboun=NULL,*ipkon=NULL,*kon=NULL;
  double *xload=NULL,*vold=NULL,*xboun=NULL,*co=NULL,*veold=NULL,*accold=NULL,
    *voldprev=NULL;
  nelemload=*nelemloadp;xload=*xloadp;sideload=*sideloadp;vold=*voldp;
  ikboun=*ikbounp;ilboun=*ilbounp;xboun=*xbounp;co=*cop;ipkon=*ipkonp;
  kon=*konp;lakon=*lakonp;veold=*veoldp;accold=*accoldp;filab=*filabp;

  if(istep==0){

    // end change DLR

   
    clock_gettime(CLOCK_MONOTONIC, &totalCalculixTimeStart);

    if(argc==1){printf("Usage: CalculiX.exe -i jobname\n");FORTRAN(stop,());}
    else{
      for(i=1;i<argc;i++){
	if(strcmp1(argv[i],"-i")==0) {
	  strcpy(jobnamec,argv[i+1]);strcpy1(jobnamef,argv[i+1],132);jin++;break;}
	if(strcmp1(argv[i],"-v")==0) {
	  printf("\nThis is Version DEVELOPMENT\n\n");
	  FORTRAN(stop,());
	}
      }
      if(jin==0){strcpy(jobnamec,argv[1]);strcpy1(jobnamef,argv[1],132);}

      /* next lines deactivated on March 18, 2020 */
    
      /*    for(i=1;i<argc;i++){
	    if(strcmp1(argv[i],"-o")==0) {
	    strcpy(output,argv[i+1]);break;}
	    }*/
    }

    putenv("CCX_JOBNAME_GETJOBNAME=jobnamec");

#ifdef BAM
    ITG lop=0,lrestart=0,kstep=1,kinc=1;
    double time[2],dtime;
    FORTRAN(uexternaldb,(&lop,&lrestart,time,&dtime,&kstep,&kinc));
#endif

    FORTRAN(openfile,(jobnamef));

    printf("\n************************************************************\n\n");
    printf("CalculiX Version DEVELOPMENT, Copyright(C) 1998-2021 Guido Dhondt\n");
    printf("CalculiX comes with ABSOLUTELY NO WARRANTY. This is free\n");
    printf("software, and you are welcome to redistribute it under\n");
    printf("certain conditions, see gpl.htm\n\n");
    printf("************************************************************\n\n");
    printf("You are using an executable made on Wed Sep 15 21:15:15 CEST 2021\n");
    fflush(stdout);

    NNEW(ipoinp,ITG,2*nentries);

    /* conservative estimate of the fields to be allocated */

    readinput(jobnamec,&inpc,&nline,&nset_,ipoinp,&inp,&ipoinpc,ithermal,&nuel_,
	      &inp_size); 

    /* saving the original input deck reading pointers */

    NNEW(ipoinp_sav,ITG,2*nentries);
    memcpy(ipoinp_sav,ipoinp,sizeof(ITG)*2*nentries);
    NNEW(inp_sav,ITG,inp_size);
    memcpy(inp_sav,inp,sizeof(ITG)*inp_size);

    /* initialization of the variables */
  
    ini_cal(jobnamec,output,fneig,kind1,kind2,itempuser,irobustdesign,&nprint,
	    neq,&mpcfree,&nbounold,&nforcold,&nloadold,&nbody_,
	    &nbodyold,&network,&nheading_,&nmpc_,&nload_,&nforc_,&nboun_,
	    &nintpoint,iperturb,&ntmat_,ithermal,&isolver,&nslavs,&nkon_,&mortar,
	    jout,&nkon,&nevtot,&ifacecount,&iplas,&npmat_,mi,&mpcend,&namtot_,
	    &iumat,&icascade,&ne1d,&ne2d,infree,&nflow,irstrt,&nener,&jrstrt,
	    &ntie_,&mcs,&nprop_,&nprop,&itpamp,&nevdamp_,&npt_,&iaxial,&inext,
	    &icontact,&nobject,&nobject_,&iit,&mpcfreeref,&isens,&namtot,&nstam,
	    &ndamp,&nef,&nk_,&ne_,&nalset_,&nmat_,&norien_,&nam_,
	    &ntrans_,&ncs_,&nstate_,&ncmat_,&memmpc_,&nprint_,energy,ctrl,alpha,
	    qaold,physcon,&istep,&istat,&iprestr,&kode,nload,&nbody,&nforc,
	    nboun,nk,&nmpc,&nam,&nzs_,nlabel,&ttime,&iheading);
  
    NNEW(set,char,81*nset_);
    NNEW(meminset,ITG,nset_);
    NNEW(rmeminset,ITG,nset_);
    NNEW(iuel,ITG,4*nuel_);

    FORTRAN(allocation,(&nload_,&nforc_,&nboun_,&nk_,&ne_,&nmpc_,&nset_,&nalset_,
			&nmat_,&ntmat_,&npmat_,&norien_,&nam_,&nprint_,mi,
			&ntrans_,set,meminset,rmeminset,&ncs_,&namtot_,&ncmat_,
			&memmpc_,&ne1d,&ne2d,&nflow,jobnamec,irstrt,ithermal,
			&nener,&nstate_,&istep,inpc,ipoinp,inp,&ntie_,&nbody_,
			&nprop_,ipoinpc,&nevdamp_,&npt_,&nslavs,&nkon_,&mcs,
			&mortar,&ifacecount,&nintpoint,infree,&nheading_,
			&nobject_,iuel,&iprestr,&nstam,&ndamp,&nef,&nbounold,
			&nforcold,&nloadold,&nbodyold,&mpcend,irobustdesign));

    SFREE(set);SFREE(meminset);SFREE(rmeminset);mt=mi[1]+1;
    NNEW(heading,char,66*nheading_);

  
    // start change DLR

  }

  // end change DLR

 
  while(istat>=0) {

    fflush(stdout);

    /* in order to reduce the number of variables to be transferred to
       the subroutines, the max. field sizes are (for most fields) copied
       into the real sizes */

    nzs[1]=nzs_;

    if((istep==0)||(irstrt[0]<0)) {
      *ne=ne_;
      nset=nset_;
      nalset=nalset_;
      nmat=nmat_;
      norien=norien_;
      ntrans=ntrans_;
      ntie=ntie_;

      /* allocating space before the first step */

      /* coordinates and topology */

      NNEW(co,double,3*nk_);
      NNEW(kon,ITG,nkon_);
      NNEW(ipkon,ITG,ne_);
      NNEW(lakon,char,8*ne_);

      /* property cards */

      if(nprop_>0){
	NNEW(ielprop,ITG,ne_);
	for(i=0;i<ne_;i++) ielprop[i]=-1;
	NNEW(prop,double,nprop_);
      }

      /* fields for 1-D and 2-D elements */

      if((ne1d!=0)||(ne2d!=0)){
	NNEW(iponor,ITG,2*nkon_);
	for(i=0;i<2*nkon_;i++) iponor[i]=-1;
	NNEW(xnor,double,36*ne1d+24*ne2d);
	NNEW(knor,ITG,24*(ne1d+ne2d)*(mi[2]+1));
	NNEW(thickn,double,2*nk_);
	NNEW(thicke,double,mi[2]*nkon_);
	NNEW(offset,double,2*ne_);
	NNEW(iponoel,ITG,nk_);
	NNEW(inoel,ITG,9*ne1d+24*ne2d);
	NNEW(rig,ITG,nk_);
	NNEW(ne2boun,ITG,2*nk_);
	if(infree[2]==0)infree[2]=1;
      }

      /* SPC's */

      NNEW(nodeboun,ITG,nboun_);
      NNEW(ndirboun,ITG,nboun_);
      NNEW(typeboun,char,nboun_+1);
      if((istep==0)||((irstrt[0]<0)&&(nam_>0)))NNEW(iamboun,ITG,nboun_);
      NNEW(xboun,double,nboun_);
      NNEW(ikboun,ITG,nboun_);
      NNEW(ilboun,ITG,nboun_);

      /* MPC's */

      NNEW(ipompc,ITG,nmpc_);
      NNEW(nodempc,ITG,3*memmpc_);
      for(i=0;i<3*memmpc_;i+=3){nodempc[i+2]=i/3+2;}
      nodempc[3*memmpc_-1]=0;
      NNEW(coefmpc,double,memmpc_);
      NNEW(labmpc,char,20*nmpc_+1);
      NNEW(ikmpc,ITG,nmpc_);
      NNEW(ilmpc,ITG,nmpc_);
      NNEW(fmpc,double,nmpc_);

      /* nodal loads */

      NNEW(nodeforc,ITG,2*nforc_);
      NNEW(ndirforc,ITG,nforc_);
      if((istep==0)||((irstrt[0]<0)&&(nam_>0)))NNEW(iamforc,ITG,nforc_);
      NNEW(idefforc,ITG,nforc_);
      NNEW(xforc,double,nforc_);
      NNEW(ikforc,ITG,nforc_);
      NNEW(ilforc,ITG,nforc_);

      /* distributed facial loads */

      NNEW(nelemload,ITG,2*nload_);
      if((istep==0)||((irstrt[0]<0)&&(nam_>0)))NNEW(iamload,ITG,2*nload_);
      NNEW(idefload,ITG,nload_);
      NNEW(sideload,char,20*nload_);
      NNEW(xload,double,2*nload_);

      /* distributed volumetric loads */

      NNEW(cbody,char,81*nbody_);
      NNEW(idefbody,ITG,nbody_);
      NNEW(ibody,ITG,3*nbody_);
      NNEW(xbody,double,7*nbody_);
      NNEW(xbodyold,double,7*nbody_);

      /* printing output */

      NNEW(prlab,char,6*nprint_);
      NNEW(prset,char,81*nprint_);

      /* set definitions */

      NNEW(set,char,81*nset);
      NNEW(istartset,ITG,nset);
      NNEW(iendset,ITG,nset);
      NNEW(ialset,ITG,nalset);

      /* (hyper)elastic constants */

      NNEW(elcon,double,(ncmat_+1)*ntmat_*nmat);
      NNEW(nelcon,ITG,2*nmat);

      /* density */

      NNEW(rhcon,double,2*ntmat_*nmat);
      NNEW(nrhcon,ITG,nmat);

      /* damping */

      if(ndamp>0){NNEW(dacon,double,nmat);}

      /* specific heat */

      NNEW(shcon,double,4*ntmat_*nmat);
      NNEW(nshcon,ITG,nmat);

      /* thermal expansion coefficients */

      NNEW(alcon,double,7*ntmat_*nmat);
      NNEW(nalcon,ITG,2*nmat);
      NNEW(alzero,double,nmat);

      /* conductivity */

      NNEW(cocon,double,7*ntmat_*nmat);
      NNEW(ncocon,ITG,2*nmat);

      /* isotropic and kinematic hardening coefficients*/

      if(npmat_>0){
	NNEW(plicon,double,(2*npmat_+1)*ntmat_*nmat);
	NNEW(nplicon,ITG,(ntmat_+1)*nmat);
	NNEW(plkcon,double,(2*npmat_+1)*ntmat_*nmat);
	NNEW(nplkcon,ITG,(ntmat_+1)*nmat);
      }

      /* linear dynamic properties */
    
      NNEW(xmodal,double,11+nevdamp_);
      xmodal[10]=nevdamp_+0.5;

      /* internal state variables (nslavs is needed for restart
	 calculations) */

      if(mortar!=1){
	NNEW(xstate,double,nstate_*mi[0]*(*ne+nslavs));
	nxstate=nstate_*mi[0]*(*ne+nslavs);
      }else if(mortar==1){
	NNEW(xstate,double,nstate_*mi[0]*(*ne+nintpoint));
	nxstate=nstate_*mi[0]*(*ne+nintpoint);
      }

      /* material orientation */

      if((istep==0)||((irstrt[0]<0)&&(norien>0))) {
	NNEW(orname,char,80*norien);
	NNEW(orab,double,7*norien);
	NNEW(ielorien,ITG,mi[2]*ne_);
      }

      /* transformations */

      if((istep==0)||((irstrt[0]<0)&&(ntrans>0))) {
	NNEW(trab,double,7*ntrans);
	NNEW(inotr,ITG,2*nk_);
      }

      /* amplitude definitions */

      if((istep==0)||((irstrt[0]<0)&&(nam_>0))) {
	NNEW(amname,char,80*nam_);
	NNEW(amta,double,2*namtot_);
	NNEW(namta,ITG,3*nam_);
      }

      if((istep==0)||((irstrt[0]<0)&&(ithermal[0]>0))) {
	NNEW(t0,double,nk_);
	NNEW(t1,double,nk_);
	if((ne1d!=0)||(ne2d!=0)){
	  NNEW(t0g,double,2*nk_);
	  NNEW(t1g,double,2*nk_);
	}
      }

      /* the number in next line is NOT 1.2357111317 -> points
	 to user input; instead it is a generic nonzero
	 initialization */

      if(istep==0){
	DMEMSET(t0,0,nk_,1.2357111319);
	DMEMSET(t1,0,nk_,1.2357111319);
      }
    
      if((istep==0)||((irstrt[0]<0)&&(ithermal[0]>0)&&(nam_>0)))NNEW(iamt1,ITG,nk_);

      if((istep==0)||((irstrt[0]<0)&&(iprestr>0)))NNEW(prestr,double,6*mi[0]*ne_);

      NNEW(vold,double,mt*nk_);
      NNEW(veold,double,mt*nk_);

      /* CFD-results */

      NNEW(vel,double,8*nef);
      NNEW(velo,double,8*nef);
      NNEW(veloo,double,8*nef);

      NNEW(ielmat,ITG,mi[2]*ne_);

      NNEW(matname,char,80*nmat);

      NNEW(filab,char,87**nlabel);

      /* tied constraints */

      if(ntie_>0){
	NNEW(tieset,char,243*ntie_);
	NNEW(tietol,double,3*ntie_);
	NNEW(cs,double,17*ntie_);
      }

      /* objectives for sensitivity analysis */

      if(nobject_>0){
	NNEW(nodedesi,ITG,nk_);
	NNEW(dgdxglob,double,nobject_*2*nk_);
	NNEW(g0,double,nobject_);
	NNEW(xdesi,double,3*nk_);
	NNEW(objectset,char,405*nobject_);
	for(i=0;i<405*nobject_;i++){objectset[i]=' ';}
      }
    
      /* temporary fields for cyclic symmetry calculations */

      if((ncs_>0)||(npt_>0)){
	if(2*npt_>24*ncs_){
	  NNEW(ics,ITG,2*npt_);
	}else{
	  NNEW(ics,ITG,24*ncs_);
	}
	if(npt_>30*ncs_){
	  NNEW(dcs,double,npt_);
	}else{
	  NNEW(dcs,double,30*ncs_);
	}
      }

      /* slave faces */

      NNEW(islavsurf,ITG,2*ifacecount+2);
    
      /* robustdesign analysis */
    
      if(irobustdesign[0]>0){
	NNEW(irandomtype,ITG,nk_);
	NNEW(randomval,double,2*nk_);
      }

    }
    else {

      /* allocating and reallocating space for subsequent steps */
      
      if((*nmethod!=4)&&(*nmethod!=5)&&(*nmethod!=8)&&(*nmethod!=9)&& 
	 ((abs(*nmethod)!=1)||(iperturb[0]<2))){
        NNEW(veold,double,mt*nk_);
      }
      else{
	RENEW(veold,double,mt*nk_);
	DMEMSET(veold,mt**nk,mt*nk_,0.);
      }
      RENEW(vold,double,mt*nk_);
      DMEMSET(vold,mt**nk,mt*nk_,0.);

      RENEW(nodeboun,ITG,nboun_);
      RENEW(ndirboun,ITG,nboun_);
      RENEW(typeboun,char,nboun_+1);
      RENEW(xboun,double,nboun_);
      RENEW(ikboun,ITG,nboun_);
      RENEW(ilboun,ITG,nboun_);

      RENEW(nodeforc,ITG,2*nforc_);
      RENEW(ndirforc,ITG,nforc_);
      NNEW(idefforc,ITG,nforc_);
      RENEW(xforc,double,nforc_);
      RENEW(ikforc,ITG,nforc_);
      RENEW(ilforc,ITG,nforc_);

      RENEW(nelemload,ITG,2*nload_);
      NNEW(idefload,ITG,nload_);
      RENEW(sideload,char,20*nload_);
      RENEW(xload,double,2*nload_);

      RENEW(cbody,char,81*nbody_);
      NNEW(idefbody,ITG,nbody_);
      RENEW(ibody,ITG,3*nbody_);
      RENEW(xbody,double,7*nbody_);
      RENEW(xbodyold,double,7*nbody_);
      for(i=7*nbodyold;i<7*nbody_;i++) xbodyold[i]=0;

      if(nam > 0) {
	RENEW(iamforc,ITG,nforc_);
	RENEW(iamload,ITG,2*nload_);
	RENEW(iamboun,ITG,nboun_);
	RENEW(amname,char,80*nam_);
	RENEW(amta,double,2*namtot_);
	RENEW(namta,ITG,3*nam_);
      }

      RENEW(ipompc,ITG,nmpc_);

      RENEW(labmpc,char,20*nmpc_+1);
      RENEW(ikmpc,ITG,nmpc_);
      RENEW(ilmpc,ITG,nmpc_);
      RENEW(fmpc,double,nmpc_);

      if(ntrans > 0){
	RENEW(inotr,ITG,2*nk_);DMEMSET(inotr,2**nk,2*nk_,0);
      }

      RENEW(co,double,3*nk_);DMEMSET(co,3**nk,3*nk_,0.);

      if(ithermal[0] != 0){
	RENEW(t0,double,nk_);DMEMSET(t0,*nk,nk_,0.);
	RENEW(t1,double,nk_);DMEMSET(t1,*nk,nk_,0.);
	if((ne1d!=0)||(ne2d!=0)){
	  RENEW(t0g,double,2*nk_);DMEMSET(t0g,2**nk,2*nk_,0.);
	  RENEW(t1g,double,2*nk_);DMEMSET(t1g,2**nk,2*nk_,0.);
	}
	if(nam > 0) {RENEW(iamt1,ITG,nk_);}
      }

    }

    /* allocation of fields in the restart file */

    if(irstrt[0]<0){
      NNEW(nodebounold,ITG,nboun_);
      NNEW(ndirbounold,ITG,nboun_);
      NNEW(xbounold,double,nboun_);
      NNEW(xforcold,double,nforc_);
      NNEW(xloadold,double,2*nload_);
      if(ithermal[0]!=0) NNEW(t1old,double,nk_); 
      NNEW(sti,double,6*mi[0]**ne);
      NNEW(eme,double,6*mi[0]**ne);
      if(nener==1)NNEW(ener,double,mi[0]**ne*2);
      if(mcs>ntie_) RENEW(cs,double,17*mcs);
      if(mortar==1){
	NNEW(pslavsurf,double,3*nintpoint);
	NNEW(clearini,double,3*9*ifacecount);
      }
			
    }

    nenerold=nener;
    nkold=*nk;

    /* opening the eigenvalue file and checking for cyclic symmetry */

    strcpy(fneig,jobnamec);
    strcat(fneig,".eig");
    cyclicsymmetry=0;
    if((f1=fopen(fneig,"rb"))!=NULL){
      if(fread(&cyclicsymmetry,sizeof(ITG),1,f1)!=1){
	printf("*ERROR reading the information whether cyclic symmetry is involved in the eigenvalue file");
	exit(0);
      }
      fclose(f1);
    }

    nmpcold=nmpc;

    /* reading the input file */


    // start change DLR

    NNEW(voldprev,double,mt*nk_);
    memcpy(voldprev,vold,sizeof(double)*mt*nk_);

    if(istep<*inputsteps){

      //   if(1){
    
      // end change DLR


      if(istep==0)mortar=-1;
      FORTRAN(calinput,(co,nk,kon,ipkon,lakon,&nkon,ne,nodeboun,ndirboun,xboun,
			nboun,ipompc,nodempc,coefmpc,&nmpc,&nmpc_,nodeforc,
			ndirforc,xforc,&nforc,&nforc_,nelemload,sideload,xload,
			nload,&nload_,&nprint,prlab,prset,&mpcfree,&nboun_,mei,
			set,istartset,iendset,ialset,&nset,&nalset,elcon,nelcon,
			rhcon,nrhcon,alcon,nalcon,alzero,t0,t1,matname,ielmat,
			orname,orab,ielorien,amname,amta,namta,&nam,nmethod,
			iamforc,iamload,iamt1,ithermal,iperturb,&istat,&istep,
			&nmat,&ntmat_,&norien,prestr,&iprestr,&isolver,fei,veold,
			timepar,xmodal,filab,jout,nlabel,&idrct,jmax,&iexpl,
			alpha,iamboun,plicon,nplicon,plkcon,nplkcon,&iplas,
			&npmat_,mi,&nk_,trab,inotr,&ntrans,ikboun,ilboun,ikmpc,
			ilmpc,ics,dcs,&ncs_,&namtot_,cs,&nstate_,&ncmat_,&iumat,
			&mcs,labmpc,iponor,xnor,knor,thickn,thicke,ikforc,ilforc,
			offset,iponoel,inoel,rig,infree,nshcon,shcon,cocon,
			ncocon,physcon,&nflow,ctrl,&maxlenmpc,&ne1d,&ne2d,&nener,
			vold,nodebounold,ndirbounold,xbounold,xforcold,xloadold,
			t1old,eme,sti,ener,xstate,jobnamec,irstrt,&ttime,qaold,
			output,typeboun,inpc,ipoinp,inp,tieset,tietol,&ntie,fmpc,
			cbody,ibody,xbody,&nbody,&nbody_,xbodyold,&nam_,ielprop,
			&nprop,&nprop_,prop,&itpamp,&iviewfile,ipoinpc,&nslavs,
			t0g,t1g,&network,&cyclicsymmetry,idefforc,idefload,
			idefbody,&mortar,&ifacecount,islavsurf,pslavsurf,
			clearini,heading,&iaxial,&nobject,objectset,&nprint_,
			iuel,&nuel_,nodempcref,coefmpcref,ikmpcref,&memmpcref_,
			&mpcfreeref,&maxlenmpcref,&memmpc_,&isens,&namtot,&nstam,
			dacon,vel,&nef,velo,veloo,ne2boun,itempuser,
			irobustdesign,irandomtype,randomval));


      // start change DLR

    }else{
      istep++;
    }

    // end change DLR

  
    SFREE(idefforc);SFREE(idefload);SFREE(idefbody);

    if(istat<0) break;
  
    /* assigning the body forces to the elements */ 

    if(nbody>0){
      ifreebody=*ne+1;
      NNEW(ipobody,ITG,2*ifreebody*nbody);
      for(k=1;k<=nbody;k++){
	FORTRAN(bodyforce,(cbody,ibody,ipobody,&nbody,set,istartset,
			   iendset,ialset,&inewton,&nset,&ifreebody,&k));
	RENEW(ipobody,ITG,2*(*ne+ifreebody));
      }
      RENEW(ipobody,ITG,2*(ifreebody-1));
    }
    
    if(irefineloop==1) {
      readnewmesh(jobnamec,nboun,nodeboun,iamboun,xboun,nload,sideload,
		  iamload,&nforc,nodeforc,iamforc,xforc,ithermal,nk,&t1,&iamt1,
		  ne,&lakon,&ipkon,&kon,istartset,iendset,ialset,set,&nset,
		  filab,&co,&ipompc,&nodempc,&coefmpc,&nmpc,&nmpc_,&labmpc,
		  &mpcfree,&memmpc_,&ikmpc,&ilmpc,&nk_,&ne_,&nkon_,&istep,
		  &nprop_,&ielprop,&ne1d,&ne2d,&iponor,&thickn,&thicke,mi,
		  &offset,&iponoel,&rig,&ne2boun,&ielorien,&inotr,&t0,&t0g,&t1g,
		  &prestr,&vold,&veold,&ielmat,irobustdesign,&irandomtype,
		  &randomval,&nalset,&nalset_,&nkon,xnor,&iaxial,
		  &network,nlabel,iuel,iperturb,&iprestr,&ntie,tieset,
		  &iparentel,ikboun,&ifreebody,&ipobody,&nbody,&iprfn,
		  &konrfn,&ratiorfn,nodempcref,coefmpcref,&memmpcref_,
		  &mpcfreeref,&maxlenmpcref,&maxlenmpc,&norien);
    }

#ifdef CALCULIX_EXTERNAL_BEHAVIOURS_SUPPORT
    for(i=0;i!=nmat;++i){
      calculix_registerExternalBehaviour(matname+80*i);
    }
#endif /* CALCULIX_EXTERNAL_BEHAVIOURS_SUPPORT */
  
    if((istep==1)&&(mortar==-1)){mortar=0;}else{icontact=1;}

    nload0=*nload;

    if(iheading==0){
      writeheading(jobnamec,heading,&nheading_);
      iheading=1;
    }

    /*    if(nheading_>=0){
      writeheading(jobnamec,heading,&nheading_);
      SFREE(heading);
      nheading_=-1;
      }*/

    if((abs(*nmethod)!=1)||(iperturb[0]<2))icascade=0;

    //    FORTRAN(writeboun,(nodeboun,ndirboun,xboun,typeboun,&nboun));

    if(istep==1) {

      SFREE(iuel);

      /* tied contact constraints: generate appropriate MPC's */

      tiedcontact(&ntie, tieset, &nset, set,istartset, iendset, ialset,
		  lakon, ipkon, kon,tietol,&nmpc, &mpcfree, &memmpc_,
		  &ipompc, &labmpc, &ikmpc, &ilmpc,&fmpc, &nodempc, &coefmpc,
		  ithermal, co, vold,&nef,&nmpc_,mi,nk,&istep,ikboun,nboun,
		  kind1,kind2);

      /* reallocating space in the first step */

      /* allocating and initializing fields pointing to the previous step */

      RENEW(vold,double,mt**nk);
      NNEW(sti,double,6*mi[0]**ne);

      /* strains */

      NNEW(eme,double,6*mi[0]**ne);

      /* residual stresses/strains */

      if(iprestr==1) {
	RENEW(prestr,double,6*mi[0]**ne);
	for(i=0;i<*ne;i++){
	  for(j=0;j<mi[0];j++){
	    for(k=0;k<6;k++){
	      sti[6*mi[0]*i+6*j+k]=prestr[6*mi[0]*i+6*j+k];
	    }
	  }
	}
      }
      else if(iprestr==2){
	RENEW(prestr,double,6*mi[0]**ne);
	for(i=0;i<*ne;i++){
	  for(j=0;j<mi[0];j++){
	    for(k=0;k<6;k++){
	      eme[6*mi[0]*i+6*j+k]=prestr[6*mi[0]*i+6*j+k];
	    }
	  }
	}
      }
      else {
	SFREE(prestr);
      }

      NNEW(nodebounold,ITG,*nboun);
      NNEW(ndirbounold,ITG,*nboun);
      NNEW(xbounold,double,*nboun);
      NNEW(xforcold,double,nforc);
      NNEW(xloadold,double,2**nload);

      /* initial temperatures: store in the "old" boundary conditions */

      if(ithermal[0]>1){
	for(i=0;i<*nboun;i++){
	  if(strcmp1(&typeboun[i],"F")==0) continue;
	  if(ndirboun[i]==0){
	    xbounold[i]=vold[mt*(nodeboun[i]-1)];
	  }
	}
      }

      /* initial temperatures: store in the "old" temperature field */

      if(ithermal[0]!=0){
	NNEW(t1old,double,*nk);
	for(i=0;i<*nk;i++) t1old[i]=t0[i];
      }

      /* element definition */

      RENEW(kon,ITG,nkon);
      RENEW(ipkon,ITG,*ne);
      RENEW(lakon,char,8**ne);

      /* property cards */

      //      if(nprop>0){
      if(nprop_>0){
	RENEW(ielprop,ITG,*ne);
	RENEW(prop,double,nprop);
      }else{
	SFREE(ielprop);SFREE(prop);
      }

      /* fields for 1-D and 2-D elements */

      if((ne1d!=0)||(ne2d!=0)){
	RENEW(iponor,ITG,2*nkon);
	RENEW(xnor,double,infree[0]);
	RENEW(knor,ITG,infree[1]);
	SFREE(thickn);
	RENEW(thicke,double,mi[2]*nkon);
	RENEW(offset,double,2**ne);
	RENEW(inoel,ITG,3*(infree[2]-1));
	RENEW(iponoel,ITG,infree[3]);
	RENEW(rig,ITG,infree[3]);
	RENEW(ne2boun,ITG,2*infree[3]);
      }

      /* set definitions */ 

      RENEW(set,char,81*nset);
      RENEW(istartset,ITG,nset);
      RENEW(iendset,ITG,nset);
      RENEW(ialset,ITG,nalset);

      /* material properties */

      RENEW(elcon,double,(ncmat_+1)*ntmat_*nmat);
      RENEW(nelcon,ITG,2*nmat);

      RENEW(rhcon,double,2*ntmat_*nmat);
      RENEW(nrhcon,ITG,nmat);

      if(ndamp>0){RENEW(dacon,double,nmat);}

      RENEW(shcon,double,4*ntmat_*nmat);
      RENEW(nshcon,ITG,nmat);

      RENEW(cocon,double,7*ntmat_*nmat);
      RENEW(ncocon,ITG,2*nmat);

      RENEW(alcon,double,7*ntmat_*nmat);
      RENEW(nalcon,ITG,2*nmat);
      RENEW(alzero,double,nmat);

      RENEW(matname,char,80*nmat);
      RENEW(ielmat,ITG,mi[2]**ne);

      /* allocating space for the state variables */

      if(mortar!=1){
	RENEW(xstate,double,nstate_*mi[0]*(*ne+nslavs));
	for(i=nxstate;i<nstate_*mi[0]*(*ne+nslavs);i++){xstate[i]=0.;}
      }else if(mortar==1){
	RENEW(xstate,double,nstate_*mi[0]*(*ne+nintpoint));
	for(i=nxstate;i<nstate_*mi[0]*(*ne+nintpoint);i++){xstate[i]=0.;}
      }

      /* next statements for plastic materials and nonlinear springs */

      if(npmat_>0){
	RENEW(plicon,double,(2*npmat_+1)*ntmat_*nmat);
	RENEW(nplicon,ITG,(ntmat_+1)*nmat);
	RENEW(plkcon,double,(2*npmat_+1)*ntmat_*nmat);
	RENEW(nplkcon,ITG,(ntmat_+1)*nmat);
      }

      /* material orientation */

      if(norien > 0) {
	RENEW(orname,char,80*norien);
	RENEW(ielorien,ITG,mi[2]**ne);
	RENEW(orab,double,7*norien);
      }
      else {
	SFREE(orname);
	SFREE(ielorien);
	SFREE(orab);
      }

      /* amplitude definitions */

      if(nam > 0) {
	RENEW(amname,char,80*nam);
	RENEW(namta,ITG,3*nam);
	RENEW(amta,double,2*namta[3*nam-2]);
      }
      else {
	SFREE(amname);
	SFREE(amta);
	SFREE(namta);
	SFREE(iamforc);
	SFREE(iamload);
	SFREE(iamboun);
      }

      if(ntrans > 0){
	RENEW(trab,double,7*ntrans);
      }
      else{SFREE(trab);SFREE(inotr);}

      //      if(ithermal[0]==0){SFREE(t0);SFREE(t1);SFREE(t0g);SFREE(t1g);}
      if(ithermal[0]==0){
	SFREE(t0);SFREE(t1);
	if((ne1d!=0)||(ne2d!=0)){
	  SFREE(t0g);SFREE(t1g);
	}
      }
      
      if((ithermal[0]==0)||(nam<=0)){SFREE(iamt1);}

      if(ncs_>0){
	RENEW(ics,ITG,ncs_);
      }else if(npt_>0){
	SFREE(ics);
      }
      //      SFREE(dcs);
      if((ncs_>0)||(npt_>0)){
	SFREE(dcs);
      }

      if(mcs>0){
	RENEW(cs,double,17*mcs);
      }else{
	SFREE(cs);
      }

      /* check for design variables */
      
      for(i=0;i<ntie;i++){
        if(strcmp1(&tieset[i*243+80],"D")==0){
          if(strcmp1(&tieset[i*243],"COORDINATE")==0){
	    icoordinate=1;
	  }
	}
      }

    }else{

      /* reallocating space in all but the first step (>1) */

      RENEW(vold,double,mt**nk);

      /* if the SPC boundary conditions were changed in the present step,
	 they have to be rematched with those in the last step. Removed SPC 
	 boundary conditions do not appear any more (this is different from
	 forces and loads, where removed forces or loads are reset to zero;
	 a removed SPC constraint does not have a numerical value any more) */
       
      NNEW(reorder,double,*nboun);
      NNEW(nreorder,ITG,*nboun);
      if(nbounold<*nboun){
	RENEW(xbounold,double,*nboun);
	RENEW(nodebounold,ITG,*nboun);
	RENEW(ndirbounold,ITG,*nboun);
      }
      FORTRAN(spcmatch,(xboun,nodeboun,ndirboun,nboun,xbounold,nodebounold,
			ndirbounold,&nbounold,ikboun,ilboun,vold,reorder,nreorder,
			mi,typeboun));
      RENEW(xbounold,double,*nboun);
      RENEW(nodebounold,ITG,*nboun);
      RENEW(ndirbounold,ITG,*nboun);
      SFREE(reorder); SFREE(nreorder);

      /* for additional forces or loads in the present step, the
	 corresponding slots in the force and load fields of the
	 previous steps are initialized */

      RENEW(xforcold,double,nforc);
      for(i=nforcold;i<nforc;i++) xforcold[i]=0;

      RENEW(xloadold,double,2**nload);
      for(i=2*nloadold;i<2**nload;i++) xloadold[i]=0;

      if(ithermal[0]!=0){
	RENEW(t1old,double,*nk);
      }

      if(nam > 0) {
	RENEW(amname,char,80*nam);
	RENEW(namta,ITG,3*nam);
	RENEW(amta,double,2*namta[3*nam-2]);
      }
    
    }

    /* reallocating fields for all steps (>=1) */

    RENEW(co,double,3**nk);

    RENEW(nodeboun,ITG,*nboun);
    RENEW(ndirboun,ITG,*nboun);
    RENEW(typeboun,char,*nboun+1);
    RENEW(xboun,double,*nboun);
    RENEW(ikboun,ITG,*nboun);
    RENEW(ilboun,ITG,*nboun);
    
    RENEW(nodeforc,ITG,2*nforc);
    RENEW(ndirforc,ITG,nforc);
    RENEW(xforc,double,nforc);
    RENEW(ikforc,ITG,nforc);
    RENEW(ilforc,ITG,nforc);

    /* temperature loading */
  
    if(ithermal[0] != 0){
      RENEW(t0,double,*nk);
      RENEW(t1,double,*nk);
      if((ne1d!=0)||(ne2d!=0)){
	RENEW(t0g,double,2**nk);
	RENEW(t1g,double,2**nk);
      }
      if(nam > 0) {RENEW(iamt1,ITG,*nk);}
    }

    RENEW(nelemload,ITG,2**nload);
    RENEW(sideload,char,20**nload);
    RENEW(xload,double,2**nload);

    RENEW(cbody,char,81*nbody);
    RENEW(ibody,ITG,3*nbody);
    RENEW(xbody,double,7*nbody);
    RENEW(xbodyold,double,7*nbody);

    RENEW(ipompc,ITG,nmpc);
    RENEW(labmpc,char,20*nmpc+1);
    RENEW(ikmpc,ITG,nmpc);
    RENEW(ilmpc,ITG,nmpc);
    RENEW(fmpc,double,nmpc);

    /* energy */

    if((nener==1)&&(nenerold==0)){
      NNEW(ener,double,mi[0]**ne*2);
      if((istep>1)&&(iperturb[0]>1)){
	printf(" *ERROR in CalculiX: in nonlinear calculations\n");
	printf("        energy output must be selected in the first step\n\n");
	FORTRAN(stop,());
      }
    }

    /* initial velocities and accelerations */

    if((*nmethod==4)||(*nmethod==5)||(*nmethod==8)||(*nmethod==9)||
       ((abs(*nmethod)==1)&&(iperturb[0]>=2))){
      RENEW(veold,double,mt**nk);
    }
    else {SFREE(veold);}

    if((*nmethod==4)&&(iperturb[0]>1)) {
      NNEW(accold,double,mt**nk);
    }

    if(nam > 0) {
      RENEW(iamforc,ITG,nforc);
      RENEW(iamload,ITG,2**nload);
      RENEW(iamboun,ITG,*nboun);
    }

    /* generate force convection elements */

    //  if(network==1){
    if(network>0){
      ne0=*ne;nkon0=nkon;nload1=*nload;
      RENEW(ipkon,ITG,*ne+*nload);
      RENEW(lakon,char,8*(*ne+*nload));
      RENEW(kon,ITG,nkon+9**nload);
      NNEW(inodesd,ITG,*nk);
      RENEW(nelemload,ITG,4**nload);
      RENEW(sideload,char,40**nload);
      
      FORTRAN(genadvecelem,(inodesd,ipkon,ne,lakon,kon,nload,
			    sideload,nelemload,&nkon,&network));
      
      SFREE(inodesd);
      RENEW(ipkon,ITG,*ne);
      RENEW(lakon,char,8**ne);
      RENEW(kon,ITG,nkon);
      RENEW(sti,double,6*mi[0]**ne);
      RENEW(eme,double,6*mi[0]**ne);
      if(iprestr>0) RENEW(prestr,double,6*mi[0]**ne);
      if(nprop>0) RENEW(ielprop,ITG,*ne);
      if((ne1d!=0)||(ne2d!=0)) RENEW(offset,double,2**ne);
      RENEW(nelemload,ITG,2**nload);
      RENEW(sideload,char,20**nload);
      RENEW(xload,double,2**nload);
      RENEW(xloadold,double,2**nload);
      if(nam>0){
	RENEW(iamload,ITG,2**nload);
	for(i=2*nload1;i<2**nload;i++)iamload[i]=0;
      }
      if(nener==1)RENEW(ener,double,mi[0]**ne*2);
      if(norien>0)RENEW(ielorien,ITG,mi[2]**ne);
      RENEW(ielmat,ITG,mi[2]**ne);
      for(i=mi[2]*ne0;i<mi[2]**ne;i++)ielmat[i]=1;
    }

    if(ntrans > 0){
      RENEW(inotr,ITG,2**nk);
    }
  
    /*   calling the user routine ufaceload (can be empty) */

    if(ithermal[1]>=2){
      NNEW(sideloadtemp,char,20**nload);
      for(i=0;i<*nload;i++){
	strcpy1(&sideloadtemp[20*i],&sideload[20*i],20);
	if((strcmp1(&sideload[20*i]," ")==0)&&
	   (strcmp1(&sideload[20*i+1]," ")!=0)){
	  strcpy1(&sideloadtemp[20*i],"F",1);
	}
      }
      FORTRAN(ufaceload,(co,ipkon,kon,lakon,nboun,nodeboun,
			 nelemload,sideloadtemp,nload,ne,nk));
      SFREE(sideloadtemp);
    }

    /* storing the undecascaded MPC's */

    if(mpcfreeref==-1){
      memmpcref_=memmpc_;mpcfreeref=mpcfree;maxlenmpcref=maxlenmpc;
      NNEW(nodempcref,ITG,3*memmpc_);
      memcpy(nodempcref,nodempc,sizeof(ITG)*3*memmpc_);
      NNEW(coefmpcref,double,memmpc_);
      memcpy(coefmpcref,coefmpc,sizeof(double)*memmpc_);
      NNEW(ikmpcref,ITG,nmpc);memcpy(ikmpcref,ikmpc,sizeof(ITG)*nmpc);
    }
 
    /* decascading MPC's only necessary if MPC's changed */

    if(((istep==1)||(ntrans>0)||(mpcend<0)||(*nk!=nkold)||(nmpc!=nmpcold))&&(icascade==0)) {
      //  if(icascade==0) {

      /* decascading the MPC's */

      printf(" Decascading the MPC's\n\n");

      callfrommain=1;
      cascade(ipompc,&coefmpc,&nodempc,&nmpc,
	      &mpcfree,nodeboun,ndirboun,nboun,ikmpc,
	      ilmpc,ikboun,ilboun,&mpcend,
	      labmpc,nk,&memmpc_,&icascade,&maxlenmpc,
	      &callfrommain,iperturb,ithermal);
    }

    /* determining the matrix structure: changes if SPC's have changed */
  
    if((icascade==0)&&(*nmethod<8)) printf(" Determining the structure of the matrix:\n");
  
    NNEW(nactdof,ITG,mt**nk);  
    NNEW(mast1,ITG,nzs[1]);
    NNEW(irow,ITG,1);
  
    if((mcs==0)||(cs[1]<0)){
      
      NNEW(icol,ITG,mt**nk);
      NNEW(jq,ITG,mt**nk+1);
      NNEW(ipointer,ITG,mt**nk);
      
      if((icascade==0)&&((*nmethod<8)||(*nmethod>10))){
	if((*nmethod==11)||(*nmethod==13)){nmethodl=2;}else{nmethodl=*nmethod;}
	mastruct(nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,nboun,ipompc,
		 nodempc,&nmpc,nactdof,icol,jq,&mast1,&irow,&isolver,neq,
		 ikmpc,ilmpc,ipointer,nzs,&nmethodl,ithermal,
		 ikboun,ilboun,iperturb,mi,&mortar,typeboun,labmpc,
		 &iit,&icascade,&network,&iexpl);
      }
      else{neq[0]=1;neq[1]=1;neq[2]=1;}
    }
    else{
      
      NNEW(icol,ITG,8**nk);
      NNEW(jq,ITG,8**nk+1);
      NNEW(ipointer,ITG,8**nk);
      
      if(*nmethod==13){nmethodl=2;}else{nmethodl=*nmethod;}
      mastructcs(nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,nboun,
		 ipompc,nodempc,&nmpc,nactdof,icol,jq,&mast1,&irow,&isolver,
		 neq,ikmpc,ilmpc,ipointer,nzs,&nmethodl,
		 ics,cs,labmpc,&mcs,mi,&mortar);
    }
  
    SFREE(ipointer);SFREE(mast1);
    if((icascade==0)&&(*nmethod<8))RENEW(irow,ITG,nzs[2]);

    /* nmethod=1: static analysis   */
    /* nmethod=2: frequency analysis  */
    /* nmethod=3: buckling analysis */
    /* nmethod=4: (linear or nonlinear) dynamic analysis */
    /* nmethod=5: steady state dynamics analysis */
    /* nmethod=6: Coriolis frequency calculation */
    /* nmethod=7: flutter frequency calculation */
    /* nmethod=8:  magnetostatics */
    /* nmethod=9:  magnetodynamics */
    /* nmethod=10: electromagnetic eigenvalue problems */
    /* nmethod=11: superelement creation */
    /* nmethod=12: sensitivity analysis  */
    /* nmethod=13: Green function calculation */
    /* nmethod=14: Robustness w.r.t. to geometric tolerances */

    if((*nmethod<=1)||(*nmethod==11)||((iperturb[0]>1)&&(*nmethod<8)))
      {
	if(iperturb[0]<2){
	
	  mpcinfo[0]=memmpc_;mpcinfo[1]=mpcfree;mpcinfo[2]=icascade;
	  mpcinfo[3]=maxlenmpc;

	  if(icascade!=0){
	    printf(" *ERROR in CalculiX: the matrix structure may");
	    printf("        change due to nonlinear equations;");
	    printf("        a purely linear calculation is not");
	    printf("        feasible; use NLGEOM on the *STEP card.");
	    FORTRAN(stop,());
	  }

	  linstatic(co,nk,&kon,&ipkon,&lakon,ne,nodeboun,ndirboun,xboun,
		    nboun,
		    ipompc,nodempc,coefmpc,labmpc,&nmpc,nodeforc,ndirforc,xforc,
		    &nforc,nelemload,sideload,xload,nload,
		    nactdof,&icol,jq,&irow,neq,&nzl,nmethod,ikmpc,
		    ilmpc,ikboun,ilboun,elcon,nelcon,rhcon,nrhcon,
		    alcon,nalcon,alzero,&ielmat,&ielorien,&norien,orab,&ntmat_,
		    t0,t1,t1old,ithermal,prestr,&iprestr,vold,iperturb,sti,nzs,
		    &kode,filab,eme,&iexpl,plicon,
		    nplicon,plkcon,nplkcon,&xstate,&npmat_,matname,
		    &isolver,mi,&ncmat_,&nstate_,cs,&mcs,&nkon,&ener,
		    xbounold,xforcold,xloadold,amname,amta,namta,
		    &nam,iamforc,iamload,iamt1,iamboun,&ttime,
		    output,set,&nset,istartset,iendset,ialset,&nprint,prlab,
		    prset,&nener,trab,inotr,&ntrans,fmpc,ipobody,ibody,xbody,
		    &nbody,
		    xbodyold,timepar,thicke,jobnamec,tieset,&ntie,&istep,&nmat,
		    ielprop,prop,typeboun,&mortar,mpcinfo,tietol,ics,
		    orname,itempuser,t0g,t1g);

	  for(i=0;i<3;i++){nzsprevstep[i]=nzs[i];}

	  memmpc_=mpcinfo[0];mpcfree=mpcinfo[1];icascade=mpcinfo[2];
	  maxlenmpc=mpcinfo[3];

	}

	else{

	  mpcinfo[0]=memmpc_;mpcinfo[1]=mpcfree;mpcinfo[2]=icascade;
	  mpcinfo[3]=maxlenmpc;

	  nonlingeo(&co,nk,&kon,&ipkon,&lakon,ne,nodeboun,ndirboun,xboun,
		    nboun,&ipompc,&nodempc,&coefmpc,&labmpc,&nmpc,nodeforc,
		    ndirforc,xforc,&nforc,&nelemload,&sideload,xload,nload,
		    nactdof,&icol,jq,&irow,neq,&nzl,nmethod,&ikmpc,&ilmpc,
		    ikboun,ilboun,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,
		    alzero,&ielmat,&ielorien,&norien,orab,&ntmat_,t0,t1,t1old,
		    ithermal,prestr,&iprestr,&vold,iperturb,sti,nzs,&kode,
		    filab,&idrct,jmax,jout,timepar,eme,xbounold,xforcold,
		    xloadold,veold,accold,amname,amta,namta,&nam,iamforc,
		    &iamload,iamt1,alpha,&iexpl,iamboun,plicon,nplicon,plkcon,
		    nplkcon,&xstate,&npmat_,&istep,&ttime,matname,qaold,mi,
		    &isolver,&ncmat_,&nstate_,&iumat,cs,&mcs,&nkon,&ener,
		    mpcinfo,output,shcon,nshcon,cocon,ncocon,physcon,&nflow,
		    ctrl,set,&nset,istartset,iendset,ialset,&nprint,prlab,
		    prset,&nener,ikforc,ilforc,trab,inotr,&ntrans,&fmpc,cbody,
		    ibody,xbody,&nbody,xbodyold,ielprop,prop,&ntie,tieset,
		    &itpamp,&iviewfile,jobnamec,tietol,&nslavs,thicke,ics,
		    &nintpoint,&mortar,&ifacecount,typeboun,&islavsurf,
		    &pslavsurf,&clearini,&nmat,xmodal,&iaxial,&inext,&nprop,
		    &network,orname,vel,&nef,velo,veloo,energy,itempuser,
		    ipobody,&inewton,t0g,t1g,&ifreebody);

	  memmpc_=mpcinfo[0];mpcfree=mpcinfo[1];icascade=mpcinfo[2];
	  maxlenmpc=mpcinfo[3];

	  for(i=0;i<3;i++){nzsprevstep[i]=nzs[i];}

	}
      }else if((*nmethod==2)||(*nmethod==13)){
      
      /* FREQUENCY ANALYSIS */
      
      if((mcs==0)||(cs[1]<0)){
#ifdef ARPACK
	  
	mpcinfo[0]=memmpc_;mpcinfo[1]=mpcfree;mpcinfo[2]=icascade;
	mpcinfo[3]=maxlenmpc;
	  
	arpack(co,nk,&kon,&ipkon,&lakon,ne,nodeboun,ndirboun,xboun,nboun,
	       ipompc,nodempc,coefmpc,labmpc,&nmpc,nodeforc,ndirforc,xforc,
	       &nforc,nelemload,sideload,xload,nload,
	       nactdof,icol,jq,&irow,neq,&nzl,nmethod,ikmpc,
	       ilmpc,ikboun,ilboun,elcon,nelcon,rhcon,nrhcon,
	       shcon,nshcon,cocon,ncocon,
	       alcon,nalcon,alzero,&ielmat,&ielorien,&norien,orab,&ntmat_,
	       t0,t1,t1old,ithermal,prestr,&iprestr,vold,iperturb,sti,nzs,
	       &kode,mei,fei,filab,
	       &iexpl,plicon,nplicon,plkcon,nplkcon,
	       &xstate,&npmat_,matname,mi,&ncmat_,&nstate_,&ener,jobnamec,
	       output,set,&nset,istartset,iendset,ialset,&nprint,prlab,
	       prset,&nener,&isolver,trab,inotr,&ntrans,&ttime,fmpc,ipobody,
	       ibody,xbody,&nbody,thicke,&nslavs,tietol,&nkon,mpcinfo,
	       &ntie,&istep,&mcs,ics,tieset,cs,&nintpoint,&mortar,&ifacecount,
	       &islavsurf,&pslavsurf,&clearini,&nmat,typeboun,ielprop,prop,
	       orname,&inewton,t0g,t1g);

	memmpc_=mpcinfo[0];mpcfree=mpcinfo[1];icascade=mpcinfo[2];
	maxlenmpc=mpcinfo[3];

	for(i=0;i<3;i++){nzsprevstep[i]=nzs[i];}

#else
	printf("*ERROR in CalculiX: the ARPACK library is not linked\n\n");
	FORTRAN(stop,());
#endif

      }else{

#ifdef ARPACK

	mpcinfo[0]=memmpc_;mpcinfo[1]=mpcfree;mpcinfo[2]=icascade;
	mpcinfo[3]=maxlenmpc;
	  
	arpackcs(co,nk,&kon,&ipkon,&lakon,ne,nodeboun,ndirboun,
		 xboun,nboun,
		 ipompc,nodempc,coefmpc,labmpc,&nmpc,nodeforc,ndirforc,xforc,
		 &nforc,nelemload,sideload,xload,nload,
		 nactdof,icol,jq,&irow,neq,&nzl,nmethod,ikmpc,
		 ilmpc,ikboun,ilboun,elcon,nelcon,rhcon,nrhcon,
		 alcon,nalcon,alzero,&ielmat,&ielorien,&norien,orab,&ntmat_,
		 t0,t1,t1old,ithermal,prestr,&iprestr,
		 vold,iperturb,sti,nzs,&kode,mei,fei,filab,
		 &iexpl,plicon,nplicon,plkcon,nplkcon,
		 &xstate,&npmat_,matname,mi,ics,cs,&mpcend,&ncmat_,
		 &nstate_,&mcs,&nkon,jobnamec,output,set,&nset,istartset,
		 iendset,ialset,&nprint,prlab,
		 prset,&nener,&isolver,trab,inotr,&ntrans,&ttime,fmpc,ipobody,
		 ibody,xbody,&nbody,&nevtot,thicke,&nslavs,tietol,mpcinfo,
		 &ntie,&istep,tieset,&nintpoint,&mortar,&ifacecount,&islavsurf,
		 &pslavsurf,&clearini,&nmat,typeboun,ielprop,prop,orname,
		 &inewton,t0g,t1g);

	memmpc_=mpcinfo[0];mpcfree=mpcinfo[1];icascade=mpcinfo[2];
	maxlenmpc=mpcinfo[3];

	for(i=0;i<3;i++){nzsprevstep[i]=nzs[i];}

#else
	printf("*ERROR in CalculiX: the ARPACK library is not linked\n\n");
	FORTRAN(stop,());
#endif

      }
    }else if(*nmethod==3){
    
#ifdef ARPACK
      arpackbu(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun,nboun,
	       ipompc,nodempc,coefmpc,labmpc,&nmpc,nodeforc,ndirforc,xforc,
	       &nforc,
	       nelemload,sideload,xload,nload,
	       nactdof,icol,jq,irow,neq,&nzl,nmethod,ikmpc,
	       ilmpc,ikboun,ilboun,elcon,nelcon,rhcon,nrhcon,
	       alcon,nalcon,alzero,ielmat,ielorien,&norien,orab,&ntmat_,
	       t0,t1,t1old,ithermal,prestr,&iprestr,
	       vold,iperturb,sti,nzs,&kode,mei,fei,filab,
	       eme,&iexpl,plicon,nplicon,plkcon,nplkcon,
	       xstate,&npmat_,matname,mi,&ncmat_,&nstate_,ener,output,
	       set,&nset,istartset,iendset,ialset,&nprint,prlab,
	       prset,&nener,&isolver,trab,inotr,&ntrans,&ttime,fmpc,ipobody,
	       ibody,xbody,&nbody,thicke,jobnamec,&nmat,ielprop,prop,
	       orname,typeboun,t0g,t1g);
#else
      printf("*ERROR in CalculiX: the ARPACK library is not linked\n\n");
      FORTRAN(stop,());
#endif
    }
    else if(*nmethod==4)
      {
	if((ne1d!=0)||(ne2d!=0)){
	  printf(" *WARNING: 1-D or 2-D elements may cause problems in modal dynamic calculations\n");
	  printf("           ensure that point loads defined in a *MODAL DYNAMIC step\n");
	  printf("           and applied to nodes belonging to 1-D or 2-D elements have been\n");
	  printf("           applied to the same nodes in the preceding FREQUENCY step with\n");
	  printf("           magnitude zero; look at example shellf.inp for a guideline.\n\n");}

	printf(" Composing the dynamic response from the eigenmodes\n\n");

	dyna(&co,nk,&kon,&ipkon,&lakon,ne,&nodeboun,&ndirboun,&xboun,nboun,
	     &ipompc,&nodempc,&coefmpc,&labmpc,&nmpc,nodeforc,ndirforc,xforc,
	     &nforc,
	     nelemload,sideload,xload,nload,
	     &nactdof,neq,&nzl,icol,irow,nmethod,&ikmpc,&ilmpc,&ikboun,&ilboun,
	     elcon,nelcon,rhcon,nrhcon,cocon,ncocon,
	     alcon,nalcon,alzero,&ielmat,&ielorien,&norien,orab,&ntmat_,&t0,
	     &t1,ithermal,prestr,&iprestr,&vold,iperturb,&sti,nzs,
	     timepar,xmodal,&veold,amname,amta,
	     namta,&nam,iamforc,iamload,&iamt1,
	     jout,&kode,filab,&eme,xforcold,xloadold,
	     &t1old,&iamboun,&xbounold,&iexpl,plicon,
	     nplicon,plkcon,nplkcon,&xstate,&npmat_,matname,
	     mi,&ncmat_,&nstate_,&ener,jobnamec,&ttime,set,&nset,
	     istartset,iendset,&ialset,&nprint,prlab,
	     prset,&nener,trab,&inotr,&ntrans,&fmpc,ipobody,ibody,xbody,&nbody,
	     xbodyold,&istep,&isolver,jq,output,&mcs,&nkon,&mpcend,ics,cs,
	     &ntie,tieset,&idrct,jmax,ctrl,&itpamp,tietol,&nalset,
	     ikforc,ilforc,thicke,&nslavs,&nmat,typeboun,ielprop,prop,orname,
	     t0g,t1g);
      }
    else if(*nmethod==5)
      {
	if((ne1d!=0)||(ne2d!=0)){
	  printf(" *WARNING: 1-D or 2-D elements may cause problems in steady state calculations\n");
	  printf("           ensure that point loads defined in a *STEADY STATE DYNAMICS step\n");
	  printf("           and applied to nodes belonging to 1-D or 2-D elements have been\n");
	  printf("           applied to the same nodes in the preceding FREQUENCY step with\n");
	  printf("           magnitude zero; look at example shellf.inp for a guideline.\n\n");}

	printf(" Composing the steady state response from the eigenmodes\n\n");

	steadystate(&co,nk,&kon,&ipkon,&lakon,ne,&nodeboun,&ndirboun,&xboun,
		    nboun,
		    &ipompc,&nodempc,&coefmpc,&labmpc,&nmpc,nodeforc,ndirforc,
		    xforc,&nforc,
		    nelemload,sideload,xload,nload,
		    &nactdof,neq,&nzl,icol,irow,nmethod,&ikmpc,&ilmpc,&ikboun,
		    &ilboun,
		    elcon,nelcon,rhcon,nrhcon,cocon,ncocon,
		    alcon,nalcon,alzero,&ielmat,&ielorien,&norien,orab,&ntmat_,
		    &t0,
		    &t1,ithermal,prestr,&iprestr,&vold,iperturb,sti,nzs,
		    timepar,xmodal,&veold,amname,amta,
		    namta,&nam,iamforc,iamload,&iamt1,
		    jout,&kode,filab,&eme,xforcold,xloadold,
		    &t1old,&iamboun,&xbounold,&iexpl,plicon,
		    nplicon,plkcon,nplkcon,xstate,&npmat_,matname,
		    mi,&ncmat_,&nstate_,&ener,jobnamec,&ttime,set,&nset,
		    istartset,iendset,ialset,&nprint,prlab,
		    prset,&nener,trab,&inotr,&ntrans,&fmpc,ipobody,ibody,xbody,
		    &nbody,
		    xbodyold,&istep,&isolver,jq,output,&mcs,&nkon,ics,cs,
		    &mpcend,
		    ctrl,ikforc,ilforc,thicke,&nmat,typeboun,ielprop,prop,
		    orname,
		    &ndamp,dacon,t0g,t1g);
      }
    else if((*nmethod==6)||(*nmethod==7))
      {

	printf(" Composing the complex eigenmodes from the real eigenmodes\n\n");

	complexfreq(&co,nk,&kon,&ipkon,&lakon,ne,&nodeboun,&ndirboun,&xboun,
		    nboun,&ipompc,&nodempc,&coefmpc,&labmpc,&nmpc,nodeforc,
		    ndirforc,xforc,&nforc,
		    nelemload,sideload,xload,nload,
		    &nactdof,neq,&nzl,icol,irow,nmethod,&ikmpc,&ilmpc,&ikboun,
		    &ilboun,
		    elcon,nelcon,rhcon,nrhcon,cocon,ncocon,
		    alcon,nalcon,alzero,&ielmat,&ielorien,&norien,orab,&ntmat_,
		    &t0,
		    &t1,ithermal,prestr,&iprestr,&vold,iperturb,&sti,nzs,
		    timepar,xmodal,&veold,amname,amta,
		    namta,&nam,iamforc,iamload,&iamt1,
		    jout,&kode,filab,&eme,xforcold,xloadold,
		    &t1old,&iamboun,&xbounold,&iexpl,plicon,
		    nplicon,plkcon,nplkcon,xstate,&npmat_,matname,
		    mi,&ncmat_,&nstate_,&ener,jobnamec,&ttime,set,&nset,
		    istartset,iendset,&ialset,&nprint,prlab,
		    prset,&nener,trab,&inotr,&ntrans,&fmpc,ipobody,ibody,xbody,
		    &nbody,
		    xbodyold,&istep,&isolver,jq,output,&mcs,&nkon,&mpcend,ics,
		    cs,
		    &ntie,tieset,&idrct,jmax,ctrl,&itpamp,tietol,&nalset,
		    ikforc,ilforc,thicke,jobnamef,mei,&nmat,ielprop,prop,orname,
		    typeboun,t0g,t1g);
      }
    else if((*nmethod>7)&&(*nmethod<12)){

      mpcinfo[0]=memmpc_;mpcinfo[1]=mpcfree;mpcinfo[2]=icascade;
      mpcinfo[3]=maxlenmpc;

      electromagnetics(&co,nk,&kon,&ipkon,&lakon,ne,nodeboun,
		       ndirboun,xboun,nboun,
		       &ipompc,&nodempc,&coefmpc,&labmpc,&nmpc,nodeforc,
		       ndirforc,xforc,
		       &nforc,&nelemload,&sideload,xload,nload,
		       nactdof,&icol,&jq,&irow,neq,&nzl,nmethod,&ikmpc,
		       &ilmpc,ikboun,ilboun,elcon,nelcon,rhcon,nrhcon,
		       alcon,nalcon,alzero,&ielmat,&ielorien,&norien,orab,
		       &ntmat_,
		       t0,t1,t1old,ithermal,prestr,&iprestr,
		       &vold,iperturb,sti,nzs,&kode,filab,&idrct,jmax,
		       jout,timepar,eme,xbounold,xforcold,xloadold,
		       veold,accold,amname,amta,namta,
		       &nam,iamforc,&iamload,iamt1,alpha,
		       &iexpl,iamboun,plicon,nplicon,plkcon,nplkcon,
		       &xstate,&npmat_,&istep,&ttime,matname,qaold,mi,
		       &isolver,&ncmat_,&nstate_,&iumat,cs,&mcs,&nkon,&ener,
		       mpcinfo,output,
		       shcon,nshcon,cocon,ncocon,physcon,&nflow,ctrl,
		       &set,&nset,&istartset,&iendset,&ialset,&nprint,prlab,
		       prset,&nener,ikforc,ilforc,trab,inotr,&ntrans,&fmpc,
		       ipobody,ibody,xbody,&nbody,xbodyold,ielprop,prop,
		       &ntie,&tieset,&itpamp,&iviewfile,jobnamec,&tietol,
		       &nslavs,thicke,
		       ics,&nalset,&nmpc_,&nmat,typeboun,&iaxial,&nload_,&nprop,
		       &network,orname,t0g,t1g);

      memmpc_=mpcinfo[0];mpcfree=mpcinfo[1];icascade=mpcinfo[2];
      maxlenmpc=mpcinfo[3];
    }

    else if(*nmethod==12){
      
      if(icoordinate==1){
        sensi_coor(co,nk,&kon,&ipkon,&lakon,ne,nodeboun,ndirboun,
		   xboun,nboun,ipompc,nodempc,coefmpc,labmpc,&nmpc,nodeforc,
		   ndirforc,xforc,&nforc,nelemload,sideload,xload,nload,
		   nactdof,icol,jq,&irow,neq,&nzl,nmethod,ikmpc,
		   ilmpc,ikboun,ilboun,elcon,nelcon,rhcon,nrhcon,
		   alcon,nalcon,alzero,&ielmat,&ielorien,&norien,orab,&ntmat_,
		   t0,t1,t1old,ithermal,prestr,&iprestr,vold,iperturb,sti,nzs,
		   &kode,filab,eme,&iexpl,plicon,
		   nplicon,plkcon,nplkcon,&xstate,&npmat_,matname,
		   &isolver,mi,&ncmat_,&nstate_,cs,&mcs,&nkon,&ener,
		   xbounold,xforcold,xloadold,amname,amta,namta,
		   &nam,iamforc,iamload,iamt1,iamboun,&ttime,
		   output,set,&nset,istartset,iendset,ialset,&nprint,prlab,
		   prset,&nener,trab,inotr,&ntrans,fmpc,ipobody,ibody,xbody,
		   &nbody,
		   xbodyold,timepar,thicke,jobnamec,tieset,&ntie,&istep,&nmat,
		   ielprop,prop,typeboun,&mortar,mpcinfo,tietol,ics,
		   &nobject,&objectset,&istat,orname,nzsprevstep,nlabel,
		   physcon,
		   jobnamef,iponor,knor,&ne2d,iponoel,inoel,&mpcend,dgdxglob,
		   g0,&nodedesi,&ndesi,&nobjectstart,&xdesi);
    
      }else{
        sensi_orien(co,nk,&kon,&ipkon,&lakon,ne,nodeboun,ndirboun,
		    xboun,nboun,ipompc,nodempc,coefmpc,labmpc,&nmpc,nodeforc,
		    ndirforc,xforc,&nforc,nelemload,sideload,xload,nload,
		    nactdof,icol,jq,&irow,neq,&nzl,nmethod,ikmpc,
		    ilmpc,ikboun,ilboun,elcon,nelcon,rhcon,nrhcon,
		    alcon,nalcon,alzero,&ielmat,&ielorien,&norien,orab,&ntmat_,
		    t0,t1,t1old,ithermal,prestr,&iprestr,vold,iperturb,sti,nzs,
		    &kode,filab,eme,&iexpl,plicon,
		    nplicon,plkcon,nplkcon,&xstate,&npmat_,matname,
		    &isolver,mi,&ncmat_,&nstate_,cs,&mcs,&nkon,&ener,
		    xbounold,xforcold,xloadold,amname,amta,namta,
		    &nam,iamforc,iamload,iamt1,iamboun,&ttime,
		    output,set,&nset,istartset,iendset,ialset,&nprint,prlab,
		    prset,&nener,trab,inotr,&ntrans,fmpc,ipobody,ibody,xbody,
		    &nbody,
		    xbodyold,timepar,thicke,jobnamec,tieset,&ntie,&istep,&nmat,
		    ielprop,prop,typeboun,&mortar,mpcinfo,tietol,ics,
		    &nobject,&objectset,&istat,orname,nzsprevstep,nlabel,
		    physcon,
		    jobnamef,iponor,knor,&ne2d,iponoel,inoel,&mpcend);
      }
    
    }

    else if(*nmethod==14){
  
      robustdesign(co,nk,&kon,&ipkon,&lakon,ne,nodeboun,ndirboun,xboun,
		   nboun,ipompc,nodempc,coefmpc,labmpc,&nmpc,nodeforc,
		   ndirforc,xforc,&nforc,nelemload,sideload,xload,nload,
		   nactdof,icol,jq,&irow,neq,&nzl,nmethod,ikmpc,ilmpc,ikboun,
		   ilboun,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
		   &ielmat,&ielorien,&norien,orab,&ntmat_,t0,t1,t1old,ithermal,
		   prestr,&iprestr,vold,iperturb,sti,nzs,&kode,filab,eme,
		   &iexpl,plicon,nplicon,plkcon,nplkcon,&xstate,&npmat_,
		   matname,&isolver,mi,&ncmat_,&nstate_,cs,&mcs,&nkon,&ener,
		   xbounold,xforcold,xloadold,amname,amta,namta,&nam,iamforc,
		   iamload,iamt1,iamboun,&ttime,output,set,&nset,istartset,
		   iendset,ialset,&nprint,prlab,prset,&nener,trab,inotr,
		   &ntrans,fmpc,ipobody,ibody,xbody,&nbody,xbodyold,timepar,
		   thicke,jobnamec,tieset,&ntie,&istep,&nmat,ielprop,prop,
		   typeboun,&mortar,mpcinfo,tietol,ics,&nobject,
		   &objectset,&istat,orname,nzsprevstep,nlabel,physcon,
		   jobnamef,iponor,knor,&ne2d,iponoel,inoel,&mpcend,
		   irobustdesign,irandomtype,randomval);
    }

    else if(*nmethod==15){
      
      crackpropagation(&ipkon,&kon,&lakon,ne,nk,jobnamec,nboun,iamboun,xboun,
		       nload,sideload,iamload,&nforc,iamforc,xforc,ithermal,t1,
		       iamt1,&co,&nkon,mi,&ielmat,matname,output,&nmat,set,
		       &nset,istartset,iendset,ialset,jmax,timepar,nelcon,
		       elcon,&ncmat_,&ntmat_,&istep,filab,nmethod);
    }

    else if(*nmethod==16){
    
      feasibledirection(&nobject,&objectset,&dgdxglob,g0,&ndesi,nodedesi,nk,
			&isolver,&ipkon,&kon,&lakon,ne,nelemload,nload,
			nodeboun,nboun,ndirboun,ithermal,co,vold,mi,&ielmat,
			ielprop,prop,&kode,nmethod,filab,&nstate_,&istep,cs,
			set,&nset,istartset,iendset,ialset,jobnamec,output,
			&ntrans,inotr,trab,orname,xdesi);         
      
    }

    SFREE(nactdof);SFREE(icol);SFREE(jq);SFREE(irow);SFREE(ipobody);

    /* check whether refinement was active */

    //          if(irefineloop==20){
    if(strcmp1(&filab[4089],"RM")==0){
      irefineloop++;

      if(irefineloop==1){
	memcpy(ipoinp,ipoinp_sav,sizeof(ITG)*2*nentries);
	memcpy(inp,inp_sav,sizeof(ITG)*inp_size);

	/* deallocating fields */

	dealloc_cal(&ncs_,&ics,&mcs,&cs,&tieset,&tietol,&co,
		    &kon,&ipkon,&lakon,&nodeboun,&ndirboun,&typeboun,&xboun,
		    &ikboun,&ilboun,&nodebounold,&ndirbounold,&xbounold,&ipompc,
		    &labmpc,&ikmpc,&ilmpc,&fmpc,&nodempc,&coefmpc,&nodempcref,
		    &coefmpcref,&ikmpcref,&nodeforc,&ndirforc,&xforc,
		    &ikforc,&ilforc,&xforcold,&nelemload,&sideload,&xload,
		    &xloadold,&cbody,&ibody,&xbody,&xbodyold,&nam,
		    &iamboun,&iamforc,&iamload,&amname,&amta,&namta,&set,
		    &istartset,&iendset,&ialset,&elcon,&nelcon,&rhcon,
		    &nrhcon,&shcon,&nshcon,&cocon,&ncocon,&alcon,
		    &nalcon,&alzero,&nprop,&ielprop,&prop,&npmat_,
		    &plicon,&nplicon,&plkcon,&nplkcon,&ndamp,&dacon,&norien,
		    &orname,&orab,&ielorien,&ntrans,&trab,&inotr,&iprestr,
		    &prestr,ithermal,&t0,&t1,&t1old,&iamt1,&ne1d,
		    &ne2d,&t0g,&t1g,irobustdesign,&irandomtype,
		    &randomval,&prlab,&prset,&filab,&xmodal,&ielmat,
		    &matname,&sti,&eme,&ener,&xstate,&vold,
		    &veold,&vel,&velo,&veloo,&iponor,&xnor,
		    &knor,&thicke,&offset,&iponoel,&inoel,&rig,
		    &ne2boun,&islavsurf,&mortar,&pslavsurf,&clearini,
		    &nobject_,&objectset,nmethod,iperturb,&irefineloop,
		    &iparentel,&iprfn,&konrfn,&ratiorfn,&heading,
		    &nodedesi,&dgdxglob,&g0,&nuel_,&xdesi);

	/* closing and reopening the output files */
	
	FORTRAN(openfile,(jobnamef));

	/* initialization of the variables */
  
	ini_cal(jobnamec,output,fneig,kind1,kind2,itempuser,irobustdesign,
		&nprint,neq,&mpcfree,&nbounold,&nforcold,&nloadold,&nbody_,
		&nbodyold,&network,&nheading_,&nmpc_,&nload_,&nforc_,&nboun_,
		&nintpoint,iperturb,&ntmat_,ithermal,&isolver,&nslavs,&nkon_,
		&mortar,jout,&nkon,&nevtot,&ifacecount,&iplas,&npmat_,mi,
		&mpcend,&namtot_,&iumat,&icascade,&ne1d,&ne2d,infree,&nflow,
		irstrt,&nener,&jrstrt,&ntie_,&mcs,&nprop_,&nprop,&itpamp,
		&nevdamp_,&npt_,&iaxial,&inext,&icontact,&nobject,&nobject_,
		&iit,&mpcfreeref,&isens,&namtot,&nstam,&ndamp,&nef,
		&nk_,&ne_,&nalset_,&nmat_,&norien_,&nam_,&ntrans_,
		&ncs_,&nstate_,&ncmat_,&memmpc_,&nprint_,energy,ctrl,alpha,
		qaold,physcon,&istep,&istat,&iprestr,&kode,nload,&nbody,&nforc,
		nboun,nk,&nmpc,&nam,&nzs_,nlabel,&ttime,&iheading);
  
	NNEW(set,char,81*nset_);
	NNEW(meminset,ITG,nset_);
	NNEW(rmeminset,ITG,nset_);
	NNEW(iuel,ITG,4*nuel_);

	FORTRAN(allocation,(&nload_,&nforc_,&nboun_,&nk_,&ne_,&nmpc_,&nset_,
			    &nalset_,&nmat_,&ntmat_,&npmat_,&norien_,&nam_,
			    &nprint_,mi,&ntrans_,set,meminset,rmeminset,&ncs_,
			    &namtot_,&ncmat_,&memmpc_,&ne1d,&ne2d,&nflow,
			    jobnamec,irstrt,ithermal,&nener,&nstate_,&istep,
			    inpc,ipoinp,inp,&ntie_,&nbody_,
			    &nprop_,ipoinpc,&nevdamp_,&npt_,&nslavs,&nkon_,&mcs,
			    &mortar,&ifacecount,&nintpoint,infree,&nheading_,
			    &nobject_,iuel,&iprestr,&nstam,&ndamp,&nef,
			    &nbounold,&nforcold,&nloadold,&nbodyold,&mpcend,
			    irobustdesign));

	SFREE(set);SFREE(meminset);SFREE(rmeminset);mt=mi[1]+1;
	NNEW(heading,char,66*nheading_);

	continue;
      }
    }

    /* reset tempuserflag */
  
    itempuser[0]=0;

    /* deleting the perturbation loads and temperatures */

    if((iperturb[0]==1)&&(*nmethod==3)) {
      nforc=0;
      *nload=0;
      nbody=0;
      if(ithermal[0]==1) {
	for(k=0;k<*nk;++k){
	  t1[k]=t0[k];
	}
      }
    }else{
      nbounold=*nboun;
      for (i=0;i<*nboun;i++) {
	nodebounold[i]=nodeboun[i];
	ndirbounold[i]=ndirboun[i];
      }
      nforcold=nforc;
      nloadold=*nload;
      nbodyold=nbody;
      
      /* resetting the amplitude to none except for time=total time amplitudes */

      if(nam > 0) {
	for (i=0;i<*nboun;i++) {
	  if(iamboun[i]>0){
	    if(namta[3*iamboun[i]-1]>0){
	      iamboun[i]=0;
	      xboun[i]=xbounold[i];}
	  }
	}
	for (i=0;i<nforc;i++){
	  if(iamforc[i]>0){
	    if(namta[3*iamforc[i]-1]>0){
	      iamforc[i]=0;
	      xforc[i]=xforcold[i];}
	  }
	}
	for (i=0;i<2**nload;i++){
	  if(iamload[i]>0){
	    if(namta[3*iamload[i]-1]>0){
	      iamload[i]=0;
	      xload[i]=xloadold[i];}
	  }
	}
	for (i=1;i<3*nbody;i=i+3){
	  if(ibody[i]>0){
	    if(namta[3*ibody[i]-1]>0){
	      ibody[i]=0;
	      xbody[7*(i-1)/3]=xbodyold[7*(i-1)/3];}
	  }
	}
	if(ithermal[0]==1) {
	  if(iamt1[i]>0){
	    if(namta[3*iamt1[i]-1]>0){
	      iamt1[i]=0;
	      t1[i]=t1old[i];}
	  }
	}
      }
    }

    /* removing the advective elements, if any */

    if(network>0){
      *ne=ne0;nkon=nkon0;
      RENEW(ipkon,ITG,*ne);
      RENEW(lakon,char,8**ne);
      RENEW(kon,ITG,nkon);
      RENEW(sti,double,6*mi[0]**ne);
      RENEW(eme,double,6*mi[0]**ne);
      if(iprestr>0) RENEW(prestr,double,6*mi[0]**ne);
      if(nprop>0) RENEW(ielprop,ITG,*ne);
      if((ne1d!=0)||(ne2d!=0)) RENEW(offset,double,2**ne);
      if(nener==1)RENEW(ener,double,mi[0]**ne*2);
      if(norien>0)RENEW(ielorien,ITG,mi[2]**ne);
      RENEW(ielmat,ITG,mi[2]**ne);

      /* reactivating the original load labels */

      for(i=*nload-1;i>=nload0;i--){
	if(strcmp2(&sideload[20*i],"                    ",20)==0){
	  iload=nelemload[2*i+1];
	  strcpy1(&sideload[20*(iload-1)],"F",1);
	}
      }

    }

    *nload=nload0;

    if((*nmethod==4)&&(iperturb[0]>1)) SFREE(accold);

    if(irstrt[0]>0){
      jrstrt++;
      if(jrstrt>=irstrt[0]){
	jrstrt=0;
	FORTRAN(restartwrite,(&istep,&nset,nload,&nforc,nboun,nk,ne,&nmpc,
			      &nalset,&nmat,&ntmat_,&npmat_,&norien,&nam,
			      &nprint,mi,&ntrans,&ncs_,&namtot,&ncmat_,&mpcend,
			      &maxlenmpc,&ne1d,&ne2d,&nflow,nlabel,&iplas,
			      &nkon,ithermal,nmethod,iperturb,&nstate_,&nener,
			      set,istartset,iendset,ialset,co,kon,ipkon,lakon,
			      nodeboun,ndirboun,iamboun,xboun,ikboun,ilboun,
			      ipompc,nodempc,coefmpc,labmpc,ikmpc,ilmpc,
			      nodeforc,ndirforc,iamforc,xforc,ikforc,ilforc,
			      nelemload,iamload,sideload,xload,elcon,nelcon,
			      rhcon,nrhcon,alcon,nalcon,alzero,plicon,nplicon,
			      plkcon,nplkcon,orname,orab,ielorien,trab,inotr,
			      amname,amta,namta,t0,t1,iamt1,veold,ielmat,
			      matname,prlab,prset,filab,vold,nodebounold,
			      ndirbounold,xbounold,xforcold,xloadold,t1old,eme,
			      iponor,xnor,knor,thicke,offset,iponoel,inoel,rig,
			      shcon,nshcon,cocon,ncocon,ics,sti,ener,xstate,
			      jobnamec,infree,prestr,&iprestr,cbody,ibody,
			      xbody,&nbody,xbodyold,&ttime,qaold,cs,&mcs,
			      output,physcon,ctrl,typeboun,fmpc,tieset,&ntie,
			      tietol,&nslavs,t0g,t1g,&nprop,ielprop,prop,
			      &mortar,&nintpoint,&ifacecount,islavsurf,
			      pslavsurf,clearini,irstrt,vel,&nef,velo,veloo,
			      ne2boun,&memmpc_,heading,&nheading_,&network));
      }
    } 


    // start change DLR

    if(ithermal[0]!=2){
      FORTRAN(calcexternalwork,(co,vold,istartset,iendset,ipkon,lakon,kon,
				ialset,&nset,set,mi,externalfaces,nelemload,
				nload,xload,sideload,delexternalwork,voldprev));
    }

    SFREE(voldprev);

    *nelemloadp=nelemload;*xloadp=xload;*sideloadp=sideload,*voldp=vold;
    *ikbounp=ikboun;*ilbounp=ilboun;*xbounp=xboun;*cop=co;*ipkonp=ipkon;
    *konp=kon;*lakonp=lakon;*veoldp=veold;*accoldp=accold;*filabp=filab;
    return;

    // end change DLR

  
  }

  FORTRAN(closefile,());

  strcpy(fneig,jobnamec);
  strcat(fneig,".frd");
  if((f1=fopen(fneig,"ab"))==NULL){
    printf("*ERROR in frd: cannot open frd file for writing...");
    exit(0);
  }
  fprintf(f1," 9999\n");
  fclose(f1);

  /* deallocating the fields
     this section is addressed immediately after leaving calinput */

  SFREE(ipoinpc);SFREE(inpc);SFREE(inp);SFREE(ipoinp);SFREE(inp_sav);
  SFREE(ipoinp_sav);

  dealloc_cal(&ncs_,&ics,&mcs,&cs,&tieset,&tietol,&co,
	      &kon,&ipkon,&lakon,&nodeboun,&ndirboun,&typeboun,&xboun,
	      &ikboun,&ilboun,&nodebounold,&ndirbounold,&xbounold,&ipompc,
	      &labmpc,&ikmpc,&ilmpc,&fmpc,&nodempc,&coefmpc,
	      &nodempcref,&coefmpcref,&ikmpcref,&nodeforc,&ndirforc,&xforc,
	      &ikforc,&ilforc,&xforcold,&nelemload,&sideload,&xload,
	      &xloadold,&cbody,&ibody,&xbody,&xbodyold,&nam,
	      &iamboun,&iamforc,&iamload,&amname,&amta,&namta,&set,
	      &istartset,&iendset,&ialset,&elcon,&nelcon,&rhcon,
	      &nrhcon,&shcon,&nshcon,&cocon,&ncocon,&alcon,
	      &nalcon,&alzero,&nprop,&ielprop,&prop,&npmat_,
	      &plicon,&nplicon,&plkcon,&nplkcon,&ndamp,&dacon,&norien,
	      &orname,&orab,&ielorien,&ntrans,&trab,&inotr,&iprestr,
	      &prestr,ithermal,&t0,&t1,&t1old,&iamt1,&ne1d,
	      &ne2d,&t0g,&t1g,irobustdesign,&irandomtype,
	      &randomval,&prlab,&prset,&filab,&xmodal,&ielmat,
	      &matname,&sti,&eme,&ener,&xstate,&vold,
	      &veold,&vel,&velo,&veloo,&iponor,&xnor,
	      &knor,&thicke,&offset,&iponoel,&inoel,&rig,
	      &ne2boun,&islavsurf,&mortar,&pslavsurf,&clearini,
	      &nobject_,&objectset,nmethod,iperturb,&irefineloop,
	      &iparentel,&iprfn,&konrfn,&ratiorfn,&heading,
	      &nodedesi,&dgdxglob,&g0,&nuel_,&xdesi);
  
#ifdef CALCULIX_MPI
  MPI_Finalize();
#endif

#ifdef CALCULIX_EXTERNAL_BEHAVIOURS_SUPPORT
  calculix_freeExternalBehaviours();
#endif /* CALCULIX_EXTERNAL_BEHAVIOURS_SUPPORT */

  clock_gettime(CLOCK_MONOTONIC, &totalCalculixTimeEnd); 

  totalCalculixTime = (totalCalculixTimeEnd.tv_sec - totalCalculixTimeStart.tv_sec) * 1e9; 
  totalCalculixTime = (totalCalculixTime + (totalCalculixTimeEnd.tv_nsec - totalCalculixTimeStart.tv_nsec)) * 1e-9;
  
  printf("________________________________________\n\n");
  
  printf("Total CalculiX Time: %lf\n", totalCalculixTime);

  printf("________________________________________\n");

  return;
      
}
